Two-point correlation functions of the diffusion-limited annihilation in one dimension 
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\ Two-point density-density correlation functions for the diffusive binary reaction system A-\- A - 

fN| . are obtained in one dimension via Monte Carlo simulation. The long-time behavior of these 

» I ' correlation functions clearly deviates from that of a recent analytical prediction of Bares and Mobilia 

CIh, [Phys. Rev. Lett. 83, 5214 (1999)]. An alternative expression for the asymptotic behavior is 

■ conjectured from numerical data. 
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Q . I. INTRODUCTION 

O , 

Nonequilibrium systems have been studied extensively during the past decades. Reaction-diffusion systems espe- 
cially have attracted much interest due to their relevance to many areas of physics, biology, economics, and so on 
^ [0,^ . A commonly used description of the dynamics of reaction-diffusion systems is the master equation. The master 
c/2 equation for the probability distribution of a many-body system is a linear equation, and can be written as an imagi- 
nary time Schrodinger equation with a (generally non-Hermitian) Hamiltonian acting on a many-particle Fock space 
The energy spectrum of the Hamiltonian contains all the information about the system. It is, in general, difficult 
to diagonalize the Hamiltonian, and many methods to treat the systems have been developed: the mean field rate 
equation, the scaling arguments [Q, and the renormalization group (RG) technique |^,^ to name only a few. 
^ I In simple situations, it is possible to find exact solutions of the model systems. A prototype of such models is the 
Q ■ one-dimensional diffusion- limited pair annihilation and/or creation of hard core particles. The model can be solved 
O I exactly only in the so-called free-fermion limit As implied in the work of Spouge JTof , the free fermion limit, for 

example, of the pair annihilation model can be reduced to a one-dimensional diffusion problem of a noninteracting 
I ' "free" particle. Any exact results for models beyond the free-fermion case would add much to our understanding of 
^ , the reaction-diffusion systems of hard core particles. 

00 ' In this context, a recent work by Bares and Mobilia is of interest. Studying the model mentioned above in the 
I general parameter space, Bares and Mobilia suggested a method to find an exact solution of a general one-dimensional 
. reaction-diffusion system applicable to the non-free-fermion case. The method relies on an analog of the Wick theorem. 
' The present authors, however, pointed out that the Wick theorem does not hold except for the free-fermion case 
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According to Ref. |12|, the method of Bares and Mobilia relies on a neglect of the higher order correlations, so the 
result of Ref. |0| can at best be regarded as a Hartree-Fock-type approximation. 

What is unknown is the validity of the approximation. In particular, the long-time behavior of the density and the 
' correlation reported in Ref. could be correct in merely a fortuitous way. To check this possibility, we performed an 
extensive Monte Carlo simulation on the system (with no bias) and compared the numerical results with the analytical 
predictions of Ref. . We find that the numerical result and the analytical prediction of Ref. for the density are 
' O ' consistent with each other in the long-time limit, while those for the two-point density-density correlation functions 
^ . clearly deviate from one another. We conjecture from the data that the two-point density-density correlation function 



Mr(t), defined below, contains an extra term with a simple analytic structure, and we estimate its coefficient. We also 
provide a heuristic argument as to why the density behavior is not modified in the Hartree-Fock-type approximation. 

II. MODEL 

In this Brief Report, we study the standard reaction-diffusion system A -\- A ^ (!> with symmetric hopping on a 
one-dimensional lattice of size L with periodic boundary conditions. Each site is either occupied by a hard core 
particle or is vacant. We restrict ourselves to the dynamic rules 

A0 with rate i, 

0A-> with rate |, (1) 
AA — > 00 with rate p. 
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where, in what follows, p is limited to be smaller than or equal to 1. To implement the above dynamics on simulation, 
we use the following algorithm: At time t, we choose one of the occupied sites randomly, and select one of the nearest 
neighbors of the chosen site with equal probability 1/2. If the selected nearest neighbor is empty, the chosen particle 
hops to that site. Otherwise, the two particles are annihilated with a probability p, and with a probability 1 — p 
nothing happens. After an attempt, the time is increased by l/N(t), where N(t) is the total particle number at a time 
t. The corresponding master equation of this algorithm can be written as = — with the Hamiltonian 



H 



E 



1 + 



1 

2^" 



(2) 



where 7 = p — 1 and ct^ represent the Pauli matrices. We regard the spin-down (-up) state as a vacuum (particle) 
state. t) is defined as \'^; t) = PiVi ^)\'^)^ where is the representation of a possible microscopic configuration, 
and P{r],t) is the probability with which the system is in state rj at time t The free-fermion limit corresponds to 
the 7 = (or p — 1) case. The above Hamiltonian, with a fully occupied initial condition, was studied in Ref. [Tl| ] and 
it was reported that the asymptotic behavior of density and correlation functions up to the subleading order become 
(we are using our notation) 



^™-4^ + {- + «'-7^}(4^' 



(3) 
(4) 



respectively, where p(t) is the density at time t, Cr(t) is the two-point (connected) correlation function, and the 
superscript BM stands for "the result of Bares and Mobilia" . 

In the following, we concentrate our attention on the two-point density-density correlation function, which is defined 

as 



Mr{t) = J ^(n„n™+^)(t) = Cr{t) + p{tf 



(5) 



where n„i = a^a:^, and (• • ■){t) represents the ensemble average at time t. The reason for studying Mr{t) rather than 
Cr{t) is twofold: First, the subleading order of the connected correlation function is the same as the leading order of 
Mr{t) [see Eqs. (^ and (|lO|)] and, in general, the study of the leading order is more reliable than that of a subleading 
order. Second, Mi{t) is an order parameter of the pair contact process with diffusion (PCPD) model ||l^, and it is 
known [Q that the PCPD model in the absorbing phase shares features with the annihilation model studied here. 
From Eqs. (^) and (||), the leading behavior of Mrit) is meant to be 



(47ri)3/2 ■ 



(6) 



To see how good Eqs. (p|) and (o) are, we present the numerical results in Sec III. 



III. NUMERICAL RESULTS 



We performed the Monte Carlo simulations for several p's (p = 1, 3/4, 1/2, and 1/4) with a system size L = 10 000. 
The initial conditions for all cases are fully occupied states as in Ref. [|ll|. For each p, we realize 10^ samples. We 
monitored the simulation up to 10^ Monte Carlo time steps. 

At first, we compare the density obtained by numerics with Eq. (|) in Fig. This fi gure clearly shows that Eq. 
(^ is in excellent agreement with the numerical results up to subleading orders. This result also implies that a system 
size of 10 000 is large enough not to show the finite size effect up to 10^ time steps. 

Next we compare the numerical results of the density-density correlation functions with Eq. (^) . For this purpose, 
in Fig. § we draw Mr{t)/M^^{t) as a function of time rather than Mr{t) itself. If Eq. (H) is correct, all data are 
expected to converge to the constant value of 1 in the long-time limit, but no such convergence is observed except for 
p = 1, i.e., the free-fermion case and for large r. Figure ^[clearly shows deviations from that anticipated by Eq. (|^). 
Hence the correct asymptotic behavior of the density-density correlation function may take the form 
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One expects A(r, p = 1) = 1 for all r. We estimate A(r, p) from the data by the least squares fits, and show them in 
Fig. js] as functions of 1/r for various p's. It appears that A(r,p) has a very simple dependence on 1/r: 

A(r,p) = ^ + 1. (8) 
r 

Interestingly, X{p) also seems to have a simple mathematical structure, i.e., \{p) oc (1 — p)/p\ see the inset of Fig. ||. 
Combining the results, we conjecture the forms of Mr{t) and C, (t) to be 

with c = 3.4 ± 0.2 in case of the fully occupied initial condition. 

In the long-time limit and at large length scales, all the p dependence is suppressed in the density and correlation 
functions. This observation is evidence of the irrelevance of 7 in the RG sense. In Sec. IV, we use this irrelevance 
argument to present a possible explanation for the partial success of Ref. [0 . 

IV. DISCUSSION AND SUMMARY 

This section discusses the success and failure of the method given in Ref. and summarizes our work. First of all, 
the incorrect result of Mr{t) is ascribed to the non-Gaussian form of the generating function (GF) [12|. The relation 
between the Gaussian form of the GF and the Wick factorization is explicitly shown in Appendix M Furthermore, 
Appendix |b| explicitly shows that Wick factorization is not always possible, even though the probability conservation 
is satisfied. From this point of view, the method of Bares and Mobilia may be regarded as an approximation scheme 
where quartic terms are neglected. As long as only the decay exponent of Mr{t) is concerned, Eq. (^) is good enough. 
The irrelevance of 7 in the RG sense makes it possible for one to obtain the correct decay exponents of density and 
correlation functions from the exactly solvable 7 = limit. 

We also resort to the irrelevance argument to explain the correct subleading behavior of the density. First note 
that the quartic coupling in the generating function, say VF, contributes to Mr for nonzero 7, but this contribution 
should be 0(t^'^/^), because W is generated by the irrelevant operator 7 and M,. ~ 0{t~^/^). Thus we expect that 
the implicit contribution of W to the density is at most 0{t~^^^) which happens to be smaller than the subleading 
order of the density decay. The Wick factorization, which neglects W, then gives a correct time dependence of the 
density up to the subleading order. We also believe that the method in Ref. would be partially successful only if 
the nonquadratic terms in the Hamiltonian are irrelevant in the RG sense. 

For completeness, we address the question of why the Wick factorization may yield the correct steady state density 
in the presence of pair creation. When the pair creation rate is not zero, the mean field rate equation for the density 
can be written as 

|=-2ey + 2e(l-pf, (11) 

with e' the annihilation rate and e the pair creation rate, using the same notation as in Ref. |ll|| . Its stationary state 
solution is Ps = 1/(1 + -^/e'/e), which interestingly is the same as the exact solution This result implies that 

the exact steady-state solution displays the characteristics of the mean field solution. Hence an exact result for the 
steady-state density from the Wick factorization is not surprising. 

In summary, we have extensively simulated the diffusion-limited annihilation model, and have shown numerically 
that the long-time behavior of the correlation functions given in Ref. ||ll[] is not correct, although the behavior of 
density is. In addition, an analytical form of the asymptotic behavior of the correlation functions is conjectured. 

ACKNOWLEDGMENTS 

This research was supported by the Catholic University of Korea research fund 2000, and by the Brain Korea 21 
Project at Seoul National University. 



3 



APPENDIX A: GENERATING FUNCTION AND WICK THEOREM 



This appendix proves that the Gaussian form of the generating function (GF) is a sufficient and necessary condition 
for Wick factorization. We restrict ourselves to the even sector of the Hilbert space as in Ref. [0|. In the following, 
we adopt the notation of Ref. ||l2|. First let us assume that the GF, discussed in Refs. and |16|] takes the Gaussian 
form 



m = exp 



^9iCg2/(9i>92) 



/ J Si3isg2 
.gi<92 



(Al) 



with the Grassmann numbers. Proving the Wick factorization for this GF is a simple combinatoric problem. One 
can easily check that the Wick theorem can be applied; hence we prove that the Gaussian form of the GF yields Wick 
factorization. 

Next we consider the following question: Does the Wick factorization imply a Gaussian form of the GF? To answer 
this, we adopt the reductio ad absurdum. Let us assume that the GF does not take a Gaussian form, but is instead 
of the form 



m = exp 



^-3i^'l2/(9l>92) 



.9l<92 



S,9l'5l32Sg3Sl?4 

9l<92<93<>34 



C91 ^92 C93 C94 (91 , 92 , , 94 ) 



(A2) 



where W (91, 927 93; 94) is the quartic coupling, and ••• stands for higher order terms. One can calculate the four 
operator correlation functions simply by differentiation, and find 



> (092 094) + (091094X092093) +1^(91,92,93,94) 



(A3) 



for (?i < 92 < 93 < 94- The presence of W prohibits the factorization. Since this result stems from the assumption of 
a non-Gaussian form of the GF, we proved that the Wick factorization implies a Gaussian form of the GF. 
Hence the Gaussian form of the GF is a sufficient and necessary condition of the Wick factorization. 



APPENDIX B: FOUR- PARTICLE ANNIHILATION MODEL 



In this appendix, we introduce a simple toy model, called the four-particle annihilation model (FPAM), defined in 
one space dimension, and analyze it to show that the probability conservation does not imply Wick factorization. The 
dynamics occurs only when four particles form a cluster. If a four-particle cluster exists, it annihilates with a rate A. 
Hopping events are not allowed. The Hamiltonian of the FPAM is 

i?FPAM = XI [^n ^n+l^«+20',;+3 " ^^n < 0",l+lO'r7+lf^rl+2^r7+2'^«+3<+3] ' (Bl) 
n 

where the periodic boundary condition is assumed. For simplicity, let us consider a fully occupied initial condition with 
a system size 4, that is, j^"; 0) — Y[i=i '^i'^l^), where |0) is the particle vacuum state and j^*; t) is the state vector as usual. 
In the momentum space |^; 0) = nq>o oJaL^|0), where possible values of g's are 7r/4 and 3tt/A. The state vector at time 
t is l^*;*) = e^"''''*!^'; 0) + (1 — e~^^*)|0). The nonvanishing two-operator correlation functions are (a-7r/407r/4) (0 ''^^'^ 
(o-37r/403Tr/4) (i) with the value e^**^* tan 7r/8 and e""''*'* cot 7r/8, respectively. The four-operator correlation function 
at time t is (a_37r/4a37r/4O-7r/4O7r/4)(0 = e"'*'*'* which is different from (a_Tr/4aTr/4)(0(o-37r/4O3^/4)(i) by e~'^^'^ — e~^^^ : 
Hence, we see a failure of the Wick theorem under the condition of probability conservation. In other words, probability 
conservation has nothing to do with Wick factorization. 
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FIG. 1. Time dependence of the density in a log-log plot. The time is measured in units of Monte Carlo steps (MCS), and 
the lattice constant is set to 1. The values of p for the solid lines are 1, 0.75, 0.5, and 0.25 from the bottom. The broken line is 
the anticipated leading behavior l/v'47rt. Inset: time dependence of p{t) — Ant for various p's (p — 0.75, 0.5, and 0.25 from 
the bottom). The broken lines are (1 — p)/{npt), and show excellent agreement with the numerical results. 



FIG. 2. Plots of t vs Mr{t) / Mf^^ it) for various p's. If Eq. (|6|) is correct, all data must converge to 1 as t goes to infinity. 
However, except for p = 1 (that is 7 = 0) , the data show a clear discrepancy. In this figure, we present the simulation results 
for r = 2 (O), 5 (A), 8 (v), and 50 (Q)- 



FIG. 3. Estimated values of A(r,p) as a function of 1/r for p = 0.25 (0)i 0.5(»), 0.75(A), and 1.0 (v)- Error bars stand for 
standard deviations of the least squares fits of Fig. H As in Fig. |l], the lattice constant is set to 1. Each data set for the same 
p seems to lie on a straight line. The fitted lines are also shown. Inset : we draw the slopes of the fitted lines as a function of 
p. These slope values also fall on a straight line (solid line), whose slope is 1.08 ± 0.03. 
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